#!/usr/bin/perl -w
##This script run BCP to get histone modification peaks. 
##need to have fragment length estimation file from HOMER


## author: Yaping Liu  lyping1986@gmail.com 
## time: 2013-9-19

#Usage:  perl wrap_BCP.pl [option] output_prefix histone.bed input.bed 

use strict;
use File::Basename;
use Getopt::Long;

sub usage {

    print STDERR "\nUsage:\n";
    print STDERR "perl wrap_BCP.pl [option] output_prefix histone.bed input.bed ... \n\n";
    
    print STDERR "  [Options]:\n\n";
	print STDERR "  --frag_len NUM : fragment length to extend the reads. transit it to BCP's -f option (Default: --frag_len -1, use HOMer to do read extension)\n\n";
	print STDERR "  --window NUM : window size to scan. transit it to BCP's -w option (Default: --window 200)\n";
    print STDERR "  --p NUM : p value. transit it to BCP's -p option (Default: --p 0.001)\n";
    exit(1);
}

my $frag_len=-1;
my $window = 200;
my $p = 0.001;
GetOptions( "frag_len=i" => \$frag_len,
			"window=i" => \$window,
			"p=f" => \$p,
);

usage() if ( scalar(@ARGV) == 0 );

my $prefix = $ARGV[0] or die "no prefix for the output";
my $histone = $ARGV[1] or die "no experimental histone modification data";
my $input = $ARGV[2] or die "no input histone modification data";

$prefix.=".BCP_HM_peaks.p-${p}.w-${window}.bed";

my $BCP = "/export/uec-gs1/laird/users/yaping/software/external/BCP_v1.1/BCP_HM";
if($frag_len == -1){
		my $homer_out = basename($histone);
		$homer_out =~ s/([\w|-]+)\S+$/$1/;
		$homer_out = dirname($histone)."/".$homer_out;
		if(!(-d "$homer_out")){
			my $cmd = "makeTagDirectory $homer_out -genome hg19 $histone\n";
			print STDERR $cmd;
			system($cmd) ==0 || die;
		}else{
			print STDERR "$homer_out exist! so existed tagInfo file will be used\n";
		}
		open(FH,"<$homer_out/tagInfo.txt") or die "no tag info file found\n";

		while(<FH>){
			if($_=~/fragmentLengthEstimate=(\d+)/){
				$frag_len=$1;
			}

		}
		close(FH);
		if($frag_len == -1){
			die "can't get fragmentLengthEstimate in $homer_out/tagInfo.txt file\n";
		}
	
}

my $cmd = "$BCP -f $frag_len -w $window -p $p -1 $histone -2 $input -3 $prefix \n";
print STDERR $cmd;
system($cmd) ==0 || die;
##BCP_HM has some bugs on chrM, so just get rid of them...
`sed -i 's/^chrM.*//g' $prefix`;

